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Abstract: Gibbs partition models arc the largest class of infinite exchange- 
able partitions of the positive integers generalizing the product form of 
the probability function of the two-parameter Poisson-Dirichlet family. Re- 
cently those models have been investigated in a Bayesian nonparametric 
approach to species sampling problems as alternatives to the Dirichlet and 
the Pitman- Yor process priors. Here we derive marginals of conditional and 
unconditional multivariate distributions arising from exchangeable Gibbs 
partitions to obtain explicit formulas for joint falling factorial moments 
of corresponding conditional and unconditional Gibbs sampling formulas. 
Our proofs rely on a known result on factorial moments of sum of non 
independent indicators. We provide an application to a Bayesian nonpara- 
metric estimation of the predictive probability to observe a species already 
observed a certain number of times. 
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1. Introduction 

Exchangeable Gibbs partitions ([16]) are the largest class of infinite exchange- 
able partitions of the positive integers generalizing the product form of ex- 
changeable partition probability function (EPPF) of the two parameter (a, 9) 
Poisson-Dirichlet partition model ([27], [31]), namely 

p a ,e{ni,...,n k ) = I |(1 - a) nj -i, (1) 

(0 + l)n-l 

for a £ (0,1), > —a, (ni,...,rih) a composition of n, 1 < k < n and 
(x) y f a = x(x + a) ■ ■ ■ {x + (y — l)a) generalized rising factorials. Their EPPF is 
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characterized by the Gibbs product form 

k 

Pa,v(ni,. ■ -,n k ) = V rhk - a) nj -i, (2) 
i=i 

for a 6 (—oo,l), and V = (V nyk ) weights satisfying the backward recursive 
relation V n . k — (n — kct)V n+ i, k + V^+i^+i, for V\\ = 1. By Theorem 12 in 
[16] each element of (2) arises as a probability mixture of extreme partitions, 
namely: Fisher's (1943) partitions ([14]) for a < 0, Ewens (0) partitions ([8], 
[22]) for a = 0, and Poisson-Kingman conditional partitions driven by the stable 
subordinator ([29]) for a G (0, 1). 

By an application of Eq. (2.6) in [30], given an infinite EPPF in the form 
(2), for each n > 1 the corresponding joint distribution of the random vector 
(iVi n • • • , ^K n ,ni K n ) of the sizes and number of the blocks in size biased order 
(i.e. in order of their least elements) is given by 

^a,v(Ni t n =nx,.. .,N Kn , n = n k ,K n = k) = (3) 

k 



V n ,kH(l-a) nr 



n k (n k + nk-i) ■ • ■ (nfe H h n±) Y[ k J= i( n j - 1) ! ' fJi 

where the combinatorial factor accounts for the number of partitions of [n] in 
which the j-th. block in order of appearance has rij elements. When the order 
of the blocks is irrelevant an alternative, more tractable coding for the joint 
distribution (3) is in exchangeable random order (cfr. Eq. (2.7) in [30]) 

iN? x = ni ,...,N e £ = n k , K n = k)= , nl ~ V n , k TT(l-a) nj _i, (4) 



a,V{ 

■lJ=l«J'- j=l 



that, from now on, we term multivariate Gibbs distribution of parameters (n, a, V). 
Corresponding Gibbs sampling formula, encoding the partition of n by the vec- 
tor of the numbers of blocks of different sizes, is obtained by the obvious change 
of variable in (4) and is given by 

P tt ,v(Ci,n = ci, . . . , C n , w = c n ) = n\V n , k TT [(1 ~," )i ~, l]C ' , (5) 

for a = X)j=i ^{ n 3 = i}, for i = 1, . . . , n, Y%=1 iCi = n anc ^ J27=i c » = Note 
that this is the general Gibbs analog of the Ewens sampling formula (cfr. [8]) 

n \0k " ^ 

^ = Cl,...,C n , n = Cn) = — ^, (6) 



encoding by the vector of counts the Dirichlet (0) partition model, ([13, 22]), 
whose EPPF is well-known to arise for a = in (1). A comprehensive reference 
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for the study of (5), also called component frequency spectrum, for general com- 
binatorial random structures is [1]. 

In this paper we study marginals of (4), both conditional and unconditional, 
in order to derive joint falling factorial moments of corresponding conditional 
and unconditional sampling formulas. Our main motivation comes from appli- 
cations in Baycsian nonparametric estimation in species sampling problems. In 
this setting, given n observations from a population of species with multiplici- 
ties of the first k species observed (ni, . . . , n^), interest may lie in conditional 
predictive estimation of quantities related to a further sample of m observations 
(cfr. e.g. [24], [25]), or in conditional estimation of some diversity index of the 
whole population, (see e.g. [5]). A common prior assumption in the Bayesian 
nonparametric approach is that the unknown relative abundances (Pi)i>\ of the 
species in the population follow a random discrete distribution belonging to the 
Gibbs family, i.e. are such that, by Kingman's correspondence (cfr. [23]), 



where (ii, . . . , %k) ranges over all ordered fc-tuples of distinct positive integers. 
This is equivalent to assume that the theoretically infinite sequence of species 
labels (Xi)i>i is exchangeable with almost surely discrete de Finetti measure 
rcprcscntable as P(-) = YaLi Pi&Yt{')> for (Pi) any rearrangement of the ranked 
frequencies (P/) satisfying (7), independent of (Yi) ~ IID H(-), for H some non 
atomic probability distribution. 

Actually the study of conditional Gibbs structures in this perspective has 
been initiated in [24] and [25] and some results for conditional falling factorial 
moments of components of (5) are in [11]. Nevertheless in those papers some 
confusion arises between conditional EPPFs, and conditional multivariate dis- 
tributions of the vector of sizes and number of the blocks in exchangeable random 
order, which heavily affects the complexity of the proofs. 

Here, after deriving marginals of conditional and unconditional multivariate 
Gibbs distributions, we obtain joint falling factorial moments of any order of 
(5), both conditional and unconditional, and explicit formulas for some distri- 
butions of interest generalizing some particular cases obtained in [11], in a direct 
way. Our analysis, besides providing a more effective technique for the study of 
Gibbs sampling formulas, with a view toward Bayesian nonparametric applica- 
tions, establishes the first systematic study of joint multivariate distributions 
arising from Gncdin-Pitman's Gibbs partition models. The paper is organized as 
follows: in Section 2 we provide marginals of (4) and, resorting to a result in [20] 
for sum of non independent indicators, derive general formulas for joint falling 
factorial moments of (5), together with some explicit marginal distributions 
and their expected values. In Section 3 we derive conditional multivariate Gibbs 
distributions and their marginals, for sizes and number of new blocks induced 
by the additional m-samplc. A complete analysis is performed for conditional 



k 



k 





3=1 



-1) 



(7) 
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Gibbs sampling formulas exploiting the same technique adopted in Section 2. In 
Section 4 we focus on multivariate Polya-like distributions arising by the con- 
ditional allocation of the additional sample in old blocks. Finally, in Section 5, 
we provide an application of marginals of multivariate Gibbs distributions to a 
Bayesian nonparametric estimation of a ra-step ahead probability to detect at 
observation n + m + 1 a species already observed a certain number of times. 

2. Marginals of multivariate Gibbs distributions 

To obtain the marginal distributions for general Multivariate Gibbs distributions 
(4) it is enough to resort to the definition of generalized central Stirling numbers 
(cfr. Eq. 1.9 and 1.19 in [30]) (see the Appendix for further details) 

i,-a nl TT ( X ~ a )n 3 -l 



S nM a ~ U\ E II 



fc! ^ 11 nA 



(m,...,nn) 3=1 



where the sum ranges over all (ni, . . . , rife) compositions of n. From now on we 
refer to (4) omitting the ex power in the notation. 

Proposition 1. Under a general Gibbs partition model (2) of parameters (a, V), 
for each n > 1 the r -dimensional marginal of (4), for < k — r < n — n j> 
is given by 

V(Ni = m,...,N r = n r ,K n = k)= (8) 



llj=l n 3- j=1 K ' ( bl ,...,b k _ r ) i=l 

for (pi, ... , bk-r) such that bi > V« and ^\ bi = n — n j ■ Multiplying and 

dividing by (n — n jY- an< ^ ~ r )' yields 

_ ~>A TT / | v Vn,k o-l-a 

~ li; L ^(" v; ; M' a)n ^ k [r] b ^u^> 

for (x) [n] = (x)(x - l)---(x-n + l). 

Corollary 1. By a known result in [16] for each model (2) the number of blocks 
K n has distribution 

P(K n = k) = V n ,kS~^ a , 

hence, conditioning (8) on K n = k yields 

P{N 1 = ni,...,N r = n r \K n = k)= (9) 



k-r 

a )bi-l 



-1,-a 



nj =1 ^!(n-E- = i^)' *[r] S nk~ a 
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independently of the specific Gibbs model. For r = k this is the general Gibbs 
analog of Eq. (41.8) in [9], and for r = 1, < fc — 1 < n — n\ and n\ 
I ..... I; ■ I 



P(iV 1=ni |X„ = fc) 



n \ (1 - S, 



»i 



l,-a, 
1 ^n-ni,fc-l 



s 



with expected value 



E(N!\K n = k) 



q— 1,— ot, — (1— a) 
" D n-l,fc-l 



/or iS n Q ' 7 generalized non-central Stirling numbers (see (57) in the Appendix). 
Marginalizing (8) with respect to K n yields 



P(iVi = ni,...,JV r = n r ) 



,k—r ' 



fc-!-=0 



2. i . Joint factorial moments of Gibbs sampling formulas 



Joint falling factorial moments for the Ewcns' sampling formula (6) of order 
(ri, . . . , r n ), for r; non negative integers and n — J2i ^ r i > 0, are in [9] (cfr. Eq. 
(41.9)) and correspond to 



E fl 



](C|,n)[ 



1=1 



nl ^)n-Y.T=i in 

(«-EF=l^)! (*)n 



n 



i=i 



Under the same conditions, the generalization to the (a, 9) Poisson-Dirichlet 
partition model (1) has been obtained is [33] and is given by 



J(Cl,n)[ 



Ll=l 



(« + «)e, 



r t — l\a 



(n-ELilny. 



n 

i=i 



[l-a) 



i-i 



n 



i 9 + a ^2 r l)n-J2ln- 
I 



In. the following Proposition we obtain the general result for the Gibbs sampling 
formula (5) by resorting to a result in [20], first established in [7] then studied 
in [21]. See also [18, 19]. 
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Proposition 2. Under a general (a, V) Gibbs partition model, joint falling 
factorial moments of the vector of counts (Ci jn , . . . , C n ^ n ) of order (r±, . . . , r n ) 
for Iri < n are given by 



E 



1=1 



t--J2i Iri 



l-l\ 



Y Vn ' kS n-Y:"lri,k-Z:, 



V ^ l ' fc-£,n=0 

(10) 

/or < /c — 5Zj =1 ri < n — Ylj—i Iri. For ri = r < |"y] and rj = /or every 
j 7^ i, t/ie r-th falling factorial moment of Ci^ n results 



n![(l - 



i—rl 



E [(Ci,n)[r]] = TmWzT^I ^.* 5 n-rl,fc-r- 



(11) 



fc-r=0 



Proof. For if„ = fc, let C; in = X^=i 1(^7 = 0- Then by a result for sum 
of non independent indicators r.v.s in Johnson & Kotz (2005, Sect. 10.2), or 
Charalambides (2005, Example 1.12), for r < n 

k 

E[(C,) M ] =E(J2HN j = l}) [r] =r\ J2 V(N ai =l,...,N ar = l), (12) 

3 = 1 (ai,...,o r ) 

where the summation is extended over all r-combinations (ai, . . . , a r ) of {1, . . . , k}. 
Since in our case the number of blocks K n is random, and the vector (Ni, . . . , N r \K n 
k) is exchangeable then, for I = 1, . . . , n, 

n — rl 



E [(Ci)[r]] = Yl E [( C l\ K » = %]] P ( A " = fc ) = 

fc-r=0 

n ~ rl fk\ 

= r\(MN 1 = l,...,N r = l\K n = k)¥(K n = k) 
Y r!f JP(JVi = /,..., N r = l,K n = k). 



k-r=0 



k-r=0 



hence 



i=i 



n-J2i Iri n 
k-J2i r,=0 2=1 llz V 



xP(JVi = l,...,JV ri =l,...,iV E;n _ rn+1 =rv..,iV E(ri = n,K n = k). (13) 

Inserting (8) in (13) the result follows. □ 

Notice that (10) generalizes the result in [11] Eq. (11), stated in terms of 
generalized factorial coefficients, (cfr. Eq. (53) and (54) in the Appendix), which 
corresponds to (11). Next Proposition generalizes Proposition 2 (Dirichlet case) 
and Proposition 4 (two parameter Poisson-Dirichlet case) in [11]. 
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Proposition 3. Under a general (a, V) Gibbs partition model, for each n > 1 
the law of C^ n , the number of blocks of size I, has distribution 

F ^i,n-x)- xl{n)x 2^ r\(l\y(n-rl-lx)l X (14) 

n — rl— xl 

Vn,kS n _' r i_ xlk _ r _ x 

-r — x—0 

for x = 0, . . . , \n/V\ , with expected value 

E(C,,„)=( (l-a^-x £ K^V-i ( 15 ) 

^ ' A;- 1=0 

and the distribution of the number of singleton species C\ >n follows from (14) 



x 

k— r— x— 



n^ n = x) = -j2 ' r '_ xV E Vn, k s-^ x<k _ r _ x . (i6) 

r=0 ' V k-r~x=0 

with expected value 

n-l 

E(C 1>n )=n E Ka^iT-r (17) 



Proof. (14) arises by the known relationship between discrete probability distri- 
butions and falling factorial moments (cfr. (58) in the Appendix). (15) follows 
from (11) for r = 1. (16) and (17) follow for / = 1 and generalize (41.10) and 
(41.11) in [9] to the entire Gibbs family. □ 



3. Conditional multivariate Gibbs distributions 



The study of conditional exchangeable random partitons, i.e. random partitions 
starting with an initial allocation of the first n natural integers in a certain num- 
ber k of blocks, has been initiated in [24] in view of proposing a Bayesian condi- 
tional nonparametric estimation of the richness of a population of species under 
priors on the unknown relative abundances belonging to the Gibbs class. (See 
also [15], Sect. 7). In [2] it has been shown that the corresponding conditional 
partition probability function, describing the conditional allocation in new and 
old blocks of integers n .+ 1, n + 2, . . . , can be obtained by a multi-step variation 
of the classical Chinese restaurant process construction (CRP) for exchangeable 
partitions, (first devised by Dubins and Pitman, see [30] Ch. 3). This varia- 
tion helps to properly place the Bayesian nonparametric approach to species 
sampling problems under Gibbs priors into the Gncdin-Pitman's exchangeable 
random partitions theoretical framework. Here we recall the multi-step CRP for 
completeness. 
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Proposition 4. (Cerquetti, 2009) Given an infinite EPPF model (2), assume 
that an unlimited numbers of groups of customers arrive sequentially in a restau- 
rant with an unlimited numbers of circular tables, each capable of sitting an 
unlimited numbers of customers. Given the placement of the first group of n 
customers in an = (rti, . . . , rij) configuration in j tables, a new group of m > 1 
customers is 

a) all seated at the j old tables in configuration m = (mi, . . . , rrij), for mi > 0, 
Ylj—i mi — m, with probability 

Pm (n) = Vr "+^ g r i( 1 - a )^+^-i = Yn±^L TT (rl . _ a)ms , (is) 

b) all seated at k new tables in configuration s = (si, . . . , Sk), for Yli=i s i = m > 
1 < k < m, Si > 1, with probability 



S( \ _ Vn+m,j+k IlLll 1 ~ ")»,-! ~ a )s,-l _ V n + m ,j+k TT,. 



k 



n+mj+k , 



V, 



i=l 



(19) 

a subset s < m of the new customers is seated at k new tables in configuration 
(si, . . . , Sfc) ttwd £/ie remaining m — s customers are seated at the old tables in 
configuration (mi, . . . , rrij) for ^Y^i=i m i = s, 1 < s < m, X)i=i s * = ,s ' 
nij > 0, Sj > 1 im£/i probability 



v n>j nLi(i-«k-i 

which, by the multiplicative property of rising factorials (4-5), simplifies to 

- Vn+m > j+k f[(n t - aU IlC 1 ~ «).«-!■ (2°) 



r, 



! = 1 



Now, as in [25], given the allocation of the first n integers in j blocks with 
multiplicities (m, ...,%), let be the number of new blocks generated by the 
additional m integers, (Si, . . . , SK m ) the vector of the sizes of the new blocks 
in exchangeable random order and S m = y Si the total number of new 
integers in new blocks. To obtain the joint conditional distribution of the vector 
(K m , S m , Si, ... , Sx m ) of the number and multiplicities of new blocks, and total 
observations in new blocks, it is enough to marginalize (20) with respect to all 
(mi, . . . , rrij) allocations of m — S m observations in old blocks, and to multiply 
for the combinatorial coefficient accounting for the number of partitions of [m] 
providing the same sizes and the same number k of new blocks and the same 
number s of integers in new blocks. We can hence state the following. 

Proposition 5. Under a general (a, V) Gibbs partition model the joint condi- 
tional distribution of (S m , K m , S\, ... , SK m ), for Si, ... , Sx m in exchangeable 
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random order, given the initial allocation of n integers in j blocks, corresponds 
to 

V(K m = k,S m =s,Si = S!,...,S Km = s k \ni,...,nj) = 



si!---s fc !fc! V„ t j \s 
or alternatively 



k 



Vn+m,j+k f m \ , ■ ■. TT/1 

'(n-ja) m _ s J[|(l 



Vn+ J l ' 3+k (n - ja) m . s f[(l a). 4 -i. (21) 



si! • • -Sklklm- s\ V n , 
Moreover, conditioning on S m , by Eq. (11) in [25], yields 

F(K m = k, Si = sx, . . . , s Km = s k \K n = j, S m = s) = (22) 



'! Vn+m,j+k 



k 



nc 1 



while conditioning on K m , by Eq. (4) in [24], eliminates the dependency on the 
specific (Ki,fc) Gibbs model as in (9) 

P(5i = si, . . . , S Km = s k \K m = k, S m = s, K n = j) = (23) 

«i nti(i 



si!---« fc !fc! ^'- Q 

Remark 1. Further results for the conditional moments of any order of if m 
and for the conditional asymptotic distribution of a proper normalization of 
K m under (a, 6) Poisson-Dirichlet partition models are in [10]. A simplified 
approach to the posterior analysis of the two-parameter model exploiting the 
deletion of classes property and the Beta-Binomial distribution of S m \K n = j is 
in [3] . A general result for conditional a diversity for Poisson-Kingman partition 
models driven by the stable subordinator ([29]) has been obtained in [4]. 
Remark 2. Notice that equations (21), (22), and (23) fix corresponding formulas 
(9), (19) and (34) in [25] which are missing the combinatorial coefficients. The 
problem in Lijoi et al. (2008) seems to follow from some confusion between con- 
ditional Gibbs EPPFs, as arising from the multistep sequential construction of 
Proposition 4, and joint conditional distributions of the corresponding random 
vectors. We stress here that an EPPF provides the probability of a particular 
partition characterized by a certain allocation in a certain number of blocks with 
certain multiplicities. This differs from the probability of the random vector of 
the multiplicities to assume that specific value, which is obtained by summing 
over all different partitions providing the same multiplicities in the same num- 
ber of blocks. The results in the following sections show that once the corrected 
formulas for the joint conditional distribution arc properly identified, the deriva- 
tion of estimators for quantities of interest in Bayesian nonparametric species 
sampling modeling simply follows by working with joint conditional marginals 
of (24). 
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The first step is to define the conditional analog of (4). 

Definition 1. Under a general (a, V) Gibbs partition model, the multivariate 
distribution of the vector {Si , . . . , Sx m j K m ), which arises by marginalizing (21) 
with respect to S m , 

Pa,v(Si =*!,..., S Km = s k ,K m = k\n) = (24) 

ml V n+m j +k ^ (n - ja) m -s rr f . v 
- Sl \... SkW . VnJ 2-i (m-s)\ 11 U a)s '- 1 

for (si,...,Sfc) : Sj £ [k,m], k £ [l,m], is termed conditional multivariate 
Gibbs distribution of parameters (m,a,j,n), for m > 1, a £ (— oo, 1) and j < ?i. 

In the next Proposition, mimicking the technique adopted in the previous 
section for the unconditional case, we derive marginals of (24) as the tools 
to obtain joint conditional falling^ factorial moments of the conditional Gibbs 
sampling formula. For Wi iTn = X^Ji = 1} this is given by the usual change 
of variable in (24), hence 

P(Wi, m = Wi, . . . , W m ,m = w m |ni,...,n J ) = (25) 
_ m!K+m,i+fc ^ (n - ja) m _ s -A- [(1 - a);-!]™ 1 



^ - tt \\ 11 



K.j ^ (m-s)! 1 \ (i\) w *Wi\ 

,J s—k 2—1 

In what follows we will resort to the convolution relation which defines non- 
central generalized Stirling numbers in terms of central generalized Stirling num- 
bers 

^^=E(jC" Q H)»- ( 2 6) 

s=k ^ ' 

see the Appendix (cfr. (55)) for further details. 

Proposition 6. Under a general (a, V) Gibbs partition model the r-dimensional 
marginal of (24), for (si, . . . , s r ) : J2i s i — s — m an d < fc — r < m — s i> 
is given by 

¥(S 1 = s 1 ,...,S r = s r ,K m = k\n)= (27) 

_ "^[IlLlC 1 ~ a )sj-l] (k - r)! Vn ±Irb j+k_ l - a -(n-j a ) 



Y[ r l = 1 s l \{m-Y J r l=1 s i )\ k\ V n j "»-EU 

Proof. Multiplying and dividing (24) by (s — Yli=i s iV- an d ( m ~ SI=i s iV- an d 
marginalizing yields 

F(6r - Sl , . . . , b r - s r , K m - k\n) - ^ ^ _ ^ ^ fc , y ^ 

Y m V- =lS " (^-ELl 5 ») ! ("-J«)m- f V- (g -SLl 3 «) ! TT/1 -A 

^ !-ok«»--)! ^ — ilw — ii (i-")^-! 

s -ELl H=k-r V ^ l=1 l > K ' (b u ...,bk- r ) LLt »=1 
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further multiplying and dividing by (k — r)\ we obtain 

TO! [IIi=l( 1 ~ (fe ~ r ) ! Vn+m,j+k x 

n- =1 s l \{m - ^21=! Si)\ k\ V nij 

m— EI=i s i 

/ TO — 



X 



^E ^r^^JCn-ia)^ 

s — ET=i Si — k — r 



and the result follows by (26). □ 

The following Proposition generalizes Theorem 2. in [11]. We adopt the no- 
tation to indicate components of (25) 

Proposition 7. Under a general (a, V) Gibbs partition model, joint falling 
factorial moments of order (n, . . . , r m ) of the conditional sampling formula (25) 
arise by an application of (27) in (13). For to — £ ; Iri > 



E 



(n) , 

(n) 

_i=l 



(28) 



m— J}; in 

E T/ c-l.-a>-(n-JQ!) 



For r/ = r < and = /or j ^ £ i/ien 



^Z.ttJHJ - / _ n| /m r y . ^ri+m,j+fcO m _ Hife _ r , 

V /• V -J n,o k _ r=0 

(29) 

which agrees with the result in Theorem 2. in [11] expressed in terms of non 
central generalized factorial numbers (cfr. (56) in the Appendix). 

Proof. By the analogy between (27) and (8) the proof moves along the same lines 
as the proof of Proposition 2, exploiting the marginals obtained in Proposition 
6. □ 

Remark 3. Notice the great computational advantage provided by the technique 
based on marginals of multivariate Gibbs distributions devised in the previous 
section with respect to the complexity of the approach adopted in [11]. (28) im- 
mediately follows as the conditional analog of the result obtained in Proposition 
2 without any need to provide a new proof. 

In [11], (cfr. Propositions 6 and 9), explicit marginals of (25) have been de- 
rived for the Dirichlet (9), the (a, 9) Poisson-Dirichlet and the Gnedin-Fisher 
(7) ([15]) partition models, . In the next Proposition wc obtain the general result 
for the entire Gibbs family thus providing the conditional analog of Proposition 
3. 
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Proposition 8. Under a general (a, V) Gibbs partition model the marginal 
distribution of (25) for x = 0, . . . , \m/l~\ corresponds to 

^^^^^ (30) 

fxl -1 / i\rlYi lr m—rl—xl 

l^- 1 -) K 1 ~ a )l-l\ v q -l,-a-(n-ja) 

X 2^ r !(mr( m _ r |-arl)! ^ ^n+^+fc^-W-*!,*-,--*- 
r=0 y ' y > k-r-x=0 

Its expected value, which provides the Bayesian nonparametric estimator under 
quadratic loss function, for the number of new species represented I times, arises 
from (29) for r = 1 

v(\y( n h - ( m \ ^ ST v , o-i--«,-(«-jo) /o-n 

H W l, m ) - I i I TT~ 2^ l7 «+™j+fc 5 m-i,fc-l ■ l^ 1 ) 

^ ' ™ J fc— 1=0 

T/ie conditional distribution of the number of new singleton species W{ rn will 
be 



V(W<- n) -<A- —— V [ ~ L > V T7 , o-l "« -("-J") 

X - Vn -3 r=0 r X >- k - r -x=0 

and a Bayesian estimator of the number of new singleton species follows from 
(31) as 



m—x , \ r m—r — x 



m — 1 

nw[%) = ^- E v n+m , j+k s-\t_~l n - ja \ 



k-l=0 



Proof. (30) arises by an application of (58), (31) follows from (29) for r = 1 and 
corresponds to Eq. (17) in [11] expressed in terms of generalized non central 
factorial numbers (cfr. (56) in the Appendix). □ 



4. Multivariate Polya-Gibbs distributions 



In this Section we focus on the conditional random allocation of the additional m 
integers in the j old blocks. First we derive the conditional joint distribution of 
the random vector (Mi tTn , . . . , M^ m , S m ) of the sizes of the m — S m observations 
falling in the j old blocks and of the total number of new observations S m falling 
in new blocks. Then, similarly to the previous sections, we move attention to the 
corresponding vector of counts and its joint falling factorial moments. From (20), 
marginalizing with respect to the partitions in new blocks, and multiplying for 
the combinatorial coefficient accounting for the number of allocations providing 
the same sizes of old blocks and the same number of total observations in new 
blocks, we obtain 

P(Afi, m = mi, . . . , M j>m = mj,S m = s\ni, ...,n 3 )= (32) 
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m - TJi \ Vn+m,j+k c -l,-a 

= u s W , S , 11^-^L— — s s, k > 

= l "*»■*• j=l fc =0 

for rrii > for i = 1, . . . ,j and X)i=i m i = m — ^m- 

Remark 4. Since the number of old blocks is fixed, (32) may be interpreted as a 
generalization of multivariate Polya distributions. If Qv n k is the conditional law, 

given (m, . . . , %•), of the vector (Pi, n , ■ ■ ■ , Pj,n,Rj,n), for P^ n = Pj\ni, . . . , n 3 
the conditional random relative abundance of the j-th species to appear, and 
Rj t n = 1 — Si=i ^*,nj then (32) turns out to be a Qy-multinomial mixture that 
we term multivariate Polya-Gibbs distribution of parameters (ni — a, . . . , rij — 
a, V). Moreover Qv will be the limit law, for m — > oo, of the random vector 

M («) q 
J "l,?ri lu j,m &m 

? • • • ) ; ; 

m mm 

where stands for a component of (32). Notice that for the two-parameter 

Poisson-Dirichlet (a, 6) model, by a result in [28], (cfr. Sect. 3.7, Corollary 20), 

(P lin ,...,P jtn ,Rj tn ) ~ Dir[ni -a,...,i%j -a,9 + ja], 

and substituting V n . k — (# + a)fc-ifa/(0 + l)n-i in (32) yields 

V a> e(M ltTn = mi, M jtm = m 3 ,S m = s|n) = 

= m! IlLiK -a) mi (9 + ja)s 
H; ; ///,:.-! (n + 9) m 

which is a proper multivariate Polya distribution of parameters (m, n\—a, . . . ,nj- 
a, + ja). 

Next Proposition provides the general marginal that we need to obtain joint 
falling factorial moments of the vector of counts corresponding to (32). 

Proposition 9. Under a general (a, V) Gibbs model, the conditional joint 
marginal distribution of the vector of the sizes {M\ t7n . . . , M rjm ) of the addi- 
tional new observations falling in the first r old blocks corresponds to 

P(Mi jTO = mi, . . . , M r>m = m r \ni, ...,nj) = (34) 

nr / \ n»-EI=i m,i 

i =1 {nj - a)„ H \ - V n+m j +k -(n-(j-r)a-£[ =1 m) 

"rc=i"»ii(m-ELi"»*) ! s F ™^' m - Er = iroi ' fe 

Proof. By (32) , the jont marginal of the first r blocks and S m is easily obtained 
as 

f(Mx, m = mi,... ,Mr, m = m r ,S m = s|n) 



nLi mi!(m - s - X)I=i mi)!s! 
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r r s y 

i=l i=l fe=0 nJ 

marginalizing with respect to S m , and multiplying and dividing by (m— Ei=i m i)' 
yields 

P(Mi )TO =m lr .., M r>m = m r \ni, ...,nj)= ^ 



s=0 



s 



which reduces to 



i=l 



and the result follows by an application of (26). □ 

Now let 0\ n ^ = J2i:n,<i H n i + Mi >m = l\ni, . . . , rij}, for I = 1, . . . , n + m, be 
the number of old blocks of size I after the allocation of the additional ?7i-samplc, 
then, to obtain the joint falling factorial moments of any order for the sampling 
formula of (32) we exploit the multivariate version of the result (12) recalled in 
the proof of Proposition 2, namely 



E 



(°£)m] = H E P ( M & =l-ni,...,M Sr = l- n r ). (35) 



For nii = I — rii, (34) specializes as 

P(Mi >m = I - m, . . . , M Ttm = I - n r \m, ...,rij)= (36) 

nl=iG - ";) ! ( m - ^ + Ei=i n ! 

v t/ c-l.-«.-(»-0'-'")«-Z}i=ini) 

fe=0 

and the one-dimensional marginal of (36) corresponds to 

, v , . m—l+rii 

ira/jir ; I \ ( m \ \ n i ~ a )l-ni T , a -l,-a,-(n-ja+a-m) 

V{M itm = i— n»|n) = I 2^ V n+m j +k S, 

t Tli J Vr 



II -J 



k=0 



(37) 

The following result easily follows from (35) as the analog of Propositions 
2 and 7. 
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Proposition 10. Under a general (a, V) Gibbs model, the joint falling factorial 
moments of the vector of the number of old blocks of different size (O^™ 
after the allocation of the additional m-sample, given the initial allocation m, . . . , 
is given by 



(") \ 
' n+m.m/? 



3 



n+m 



(II( £)[n]) 



Z=l 



n+m 

n 



E 



1=1 (S ri ,...,S r „ 



m-J2i in+Ei Elii "s, T/ 

'n+mj+t „-l,-tt,-(>i-(j-Ein)a-£i Ei"€<) 

^ K, "»-Ej'n+EjErii"«i.* 

/or E n = (£i,...,£ ri ),..., 3 r „ +w> = (fe,r|-r„ +m ,---)^r = + rJ' & 1 - 

and each 3 n ranges over all the combinations of n elements of j . For n = r 
and Tj = for j I, then 



E 



(o£) H 



E 



(fi, -,£,-) 



n«=i ( z - n^)!(m -rZ + £ i=1 n «> 



E 



Vn+m,j+k q—h—ot, — (n-(J-r)a—J2i=i "«i) 
T/ m— (r+Ei-i "e- > fc 



(38) 



for £j : < I, which agrees with the result in Theorem 1. in [11]. 

Next Proposition generalizes the results in Proposition 5 (two-parameter 
Poisson-Dirichlet case) and Proposition 9 (one parameter Gncdin-Fishcr case 
[15]) in [11] to the entire (a, V) Gibbs family. 

Proposition 11. Under a general (a, V) Gibbs model, from (38) and (58), the 
conditional marginal law of 0\zl is given by 



E 



I 

E 

m! 



1-2/ 



(-ir(r + y)!_L_ 



y\r\ 



r+y 



t) m=f a - ^)!(m - rl - t/Z + n,J! f = * 



E 

fc=0 



5" 



-1 -a,-(n-(j-i'-9)a-Efc 1 1 ' 
i— Ir— ij/+ET=i "Si 



A. Cerquetti/ Marginals of multivariate Gibbs distributions 



16 



Its expected value, which plays the role of the Bayesian nonparametric estimator, 
under quadratic loss function, of the number of old species represented I times, 
follows from (37) as 



E(Og) = E ( Y, ifa + M i.m = l\ni, •..,%)) = (39) 
= E(l(Mi, m = l-n i \m,...,n j ))= ^ P(M i)m = I -m\ni,...,nj) = 

i:ni<l i:ni<l 

E, \ I \ m—l+rii 

I m \ [rii — ct)i_ ni y-^ i j _ Q ,._(„_ :)Q , +Q ,_ n .) 
U-nJ K~ ^ v n+m , j+k t, m _ l+ni<k 
i:m<l v l/ nj fc=o 

or from (38) /or r = 1 and agrees with Eq. (15) in [ll]. 

Remark 5. Relying on the technique presented in this paper, falling factorial 

(n) 

r-th moments of Z\ ' , the total number of old and new blocks of size I after the 
allocation of the additional m-sample, as derived in Th. 3 in [11] by means of a 
very complex procedure, may be obtained in a straightforward way by the full 
conditional joint distribution 

f(Si =s 1 ,...,S k = s k , S m = s, K m = k, M 1>m = mi, . . . , M jiTn = m^n) = 

I i/ k j 

ml 



i=l 



Multiplying for the way to choose t blocks among the old and r — t among the 
new for every t, from (29) and (38) we get 



E 



£fl> £ 



(l-a) z _i] r 'n i=1 K, -a)i-n £i 



v {u Xi h) nuo - n 6 )!(/!)'-'(m - « + e; =1 - (»• - *)*)' 

E*n+m,fc+j ^— 1,— a, — (n— ia- 5^ n {j+* a ) 
~Y . m—rl+^ J n^ i ,k — r+t 

k-r+t=0 n ' J 

which agrees with Theorem 3. in [11]. 

In the next section we provide one more example of the importance of working 
with marginals of conditional multivariate Gibbs distributions in the implemen- 
tation of the Bayesian nonparametric approach to species sampling problems 
under Gibbs priors. 
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5. Bayesian nonparametric estimation of the probability to observe 
a species of a certain size 

In species sampling problems, particularly in ecology or genomics, given a basic 
sample (jii, . . . , nj), interest may be in estimating the probability to observe at 
step n + m + 1 a species already represented I times both belonging to an old 
species or to a new species eventually arising in the m-additional sample which 
is still not observed. This is the topic of a recent paper by Favaro et al. (2012b) 
and can be seen as a generalization of the problem of estimating the discovery 
probability, i.e. the probability to discover a new species, not represented in 
the previous n + m observations. A Bayesian nonparametric estimator of the 
discovery probability under general (a, V) Gibbs partition models has been first 
derived in [24]. 

In this Section we show that working with marginals of conditional Gibbs mul- 
tivariate distributions greatly simplifies the derivation of the results obtained 
in [12], thus providing another example of the importance of the technique pro- 
posed in this paper. 

First recall that by sequential construction of exchangeable partitions, the 
probability to observe an old species observed I times in the basic n-sample at 
observation n + 1, easily follows by one-step prediction rules for general Gibbs 
EPPFs (sec e.g. [30]). For cj,„ = Y% =1 l{n t = I}, for / = 1, . . . , n then 

/ \ p(ni,...,Z + l,...,nj) Vn+i,j n n 

Pl,n{ni, ...,Tlj) = Cz,„ : — = Q,„— (f - a). 

p{ni,...,l,...,nj) Kj 

Given a basic sample (rai, . . . ,rij), but assuming as in [12] an intermediate 
m-sample still to be observed, the probability to observe a species represented 
I times among new species at observation n + m + 1 will be a random variable, 
namely 

v) = Y^±m±^ {1 _ a)w M, (40) 

for K m the random number of new species induced by the additional sample 
and the random number of species represented / times in the additional 

sample given the basic sample. 

In the following Proposition we show how the Bayesian nonparametric esti- 
mator, under quadratic loss function, of (40), (cfr. Theorem 2. in [12]), may be 
obtained in few elegant steps. 

Proposition 12. Under a general (a, V) Gibbs partition model, for = 

Sfca = l\Kn = j} the Bayesian nonparametric estimator of P™ e ~^ 1 n I +1 (a, V) 
is given by 

(Si,...,S Km ,K m \K n =i) ~Tr (I ~ a)Wi ( J = 

\ Vn+m,j+K m J 
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( l ~ a > 1^ — l l )0--u)l-iS m \ k _i '■ (41) 



fc-l=0 



""J 



Proof. Let f(K m ) = v^) 3 + + K Km then, by definition of W^, 

^(Si s Km ,K m \K n =j) -T7 (l-aWJ = 

\ Vn+mJ+K m J 
m-l + 1 / k \ 



fe=l 



m—l 



xF{K m = k\K n = j) = 



= (l-a) J2 f(k)kP(Si = l\K m = k,K n =j)F(K m = k\K n =j). (42) 

fe-l=0 

Now, specializing (27), 

P(5i =/,..., S r = i, A™ = fc| K n = j) = 
ml (k- r)\ [(1 - a)t-i] r K+mj+fc ^_i,_a,_( n _j Q ) 



m-rZ! (/!) r fe! V n . 3 "m-ri,k- r 

and inserting the marginal for r — 1 in (42), the result follows. □ 

By analogous approach we provide a straightforward derivation for the Baycsian 
nonparametric estimator for the probability to observe a species represented / 
times among the old species, namely 

v) = Y^n±M±^ {l _ o)0 <w 

yn-\-m,j-\-K m 

Proposition 13. Under a general (a, V) Gibbs partition model, for = 

X)i=i + A^i.m = i/ien a Bayesian nonparametric estimator 

under quadratic loss function of P"iJ~™ +1 (a,V n j) is given by 

E (Ml , m ,... MjMni ,..., nj) - a)0W) = (43) 

\ Vn,m,j+K m J 



n \P ( m \r<- \ Vn+m+l,j+k a -l,-a,-(n-ja+S-a) 

= {l-a)}^md ,_J(£-a)i-e 2^ tt" ^m-i+s,fc 

e=i \ fc=o " J 

Proo/. Let f(K m ) = ^±1^ , then 

E(M 1 . m ,...,M 3 -. m ,K m |n 1 ,...,n 3 -) "fT (* ~ a )°I,m J = 
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m 

= G-a)X)/(fc): 



k=0 

j 



X ^(M 1 , m ,...,M j , m \K m =k,n u ...,n j ) 1 ^ + Af»,m = i|-Km = fe> «1 5 ■ • ■ 

\t=l > 

xV(K m = k\K n =j), 

which is equal to 

m 

= i l - a )^2 f( k ) X] P ( Af i.™ = ' - "ii-^m = k \ n l, Tlj) 
k—0 i:rti<l 

and by (37) 

i'; 1 \ " I m \t \ \ ^ fn+m+lj+fe c 

= (!-<*) 2^ 2^ — tT- — 5 



»ri+m+lji+fe 1 ,-a.-(n-ja-Tij+Q:) 



and the result follows. □ 



Appendix A 

This Appendix contains some basic facts on rising and falling factorial numbers, 
partitions and compositions of the natural integers, together with known results 
and definitions of generalized central and non central Stirling numbers that are 
exploited in the proofs and derivations all over the paper. The main reference is 
[30]. Additionally, to facilitate the reading of the results contained in [24, 25, 11] 
and [12], the relationship between central and non central generalized factorial 
coefficients and generalized Stirling numbers is reported. 



A . 1 . Generalized rising factorials 

For n = 0,1,2,..., and arbitrary real x and h, {x) n ^h denotes the nth factorial 
power of x with increment h (also called generalized rising factorial) 

n-l 

(x)„ tA := x(x + h)--- (x + (n - l)h) = J[ {x + ih) = h n (x/h) n , (44) 

i=0 

where (x) n stands for (x) n ^i, and (x)^ = x h , for which the following multi- 
plicative law holds 

{x) n +rfh = {x) n fh(x + nh) r ^ h . (45) 

From e.g. [26] (cfr. eq. 2.41 and 2.45) a binomial formula also holds, namely 

(x + y) n -th = ( D ( x )kth(v)n-kth, (46) 



fc=0 
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as well as a generalized version of the multinomial theorem, i.e. 
p l p 

d>W = E ni i..'. n (47) 

For nij > 0, for every j, and £^ = m, an application of (45) yields 

{Zj) nj+mj - x = (zj^-^Zj + rrij - (48) 

and by (47) 

„! p p p 

/fcl . ' ' * / l/<ri • 

p p 

= n^)m^-x("i + E«,- -p) n , 
3=1 3=1 

which makes unnecessary the proof of Lemma 1 in [25] . 
j4.2. Partitions and compositions 

A partition of the finite set [n] = (1, . . . , n) into k blocks is an unordered collec- 
tion of non-empty disjoint sets {A\, . . . , A^} whose union is [n], where the blocks 
Ai are assumed to be listed in order of appearance, i.e. in the order of their least 
elements. The sequence (\A\\, . . . , \Ak\) of the sizes of blocks, (m, . . . ,Hfe), de- 
fines a composition of n, i.e. a sequence of positive integers with sum n and 
■p£, denotes the space of all partitions of [n] with k blocks. From [30] (cfr. cq. 
(1.9)) the number of ways to partition [n] into k blocks and assign each block a 
W combinatorial structure such that the number of W-structures on a set of j 
elements is Wj , in terms of sum over compositions of n into k parts is given by 

^(«.) = g E fl^T> (^) 

(ni,...,nfc) 2—1 

where B n ^(w,) is a polynomial in variables w±, . . . , w n -k+i known as the (n, fc)th 
partial Bell polynomial. 

A. 3. Generalized Stirling numbers 

(For a comprehensive treatment see [17], see also [30] Ex. 1.2.7). For arbitrary 
distinct reals r\ and (3, these are the connection coefficients S 1 ^ defined by 

n 

(x)n±r, = E S n,ki X )U0 ( 50 ) 
fe=0 
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and correspond to 

SZi = Bn,k((J3-Ti).-nn)> 
where (x) n \,h are generalized falling factorials and (x) n i-h = (a;)ntfr> while 
{x)n\i = { x )[n}- Hence for r\ = — 1, /3 = —a, and a € (— oo, 1), 5~ fc ' _Q is 
defined by 

n 

=53 S, n,fe _a ( a; )fet«> ( 51 ) 

and for w ru = (1 — a) ni _i and a G [0, 1), equation (49) yields 

k 

B n , k ((i - o)._i) = ]r n^ 1 = 

(ni,...,n fc ) i=l 

In [24, 25, 11, 12] the treatment is in term of generalized factorial coefficients, 
which are the connection coefficients C" k defined by 



(ay) n = Y^Cn,k(yh, (53) 

fc=0 

(cfr. [6]). From (44) and (51), if x = ya then 

n n 
k=0 k=0 

hence 

C fe '" Q = (54) 
The representation (37) in [25], ([32]), also holds for generalized Stirling numbers 
with the obvious changes (cfr. e.g. [30], Eq. 3.19). Additionally, specializing 
formula (16) in [17], the following convolution relation holds, which defines non- 
central generalized Stirling numbers 

C = E(j^(-7U ( 55 ) 

s=k ^ ' 

and by (54), 



C2 = a k S-^ = J2 [ s )C?A-l)n- s . (56) 

Hence the following variation of equation (38) in [25] defines non-central gener- 
alized Stirling numbers as the connection coefficients S~ k ' ~ a ' 1 such that 

n n 

{ya - i) n = S~^~ an a k {y) k = ^ S~\ ~ an {ya) k ^ a . (57) 
fe=o fe=o 
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A. 4- Factorial moments and discrete distributions 

(Cfr. e.g. [20]). Falling factorial moments of a discrete r.v. X provide the distri- 
bution function by the relationship 

F(X = x) = J2^-n(X) [x+r] ] (58) 

r>0 

and standard moments by the definition as connection coefficients of the Stirling 
numbers of the second kind S r 'j , (cfr. Eq. (50)) 

r 

E[(XY}=Y l S° r ;jE[(X) [r] ], (59) 

3=0 
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